Data analysis workflow examples

library(OmniOmics)
library(MSnbase)
library(ggplot2)
library(plotly)
library(xcms)
library(CAMERA)
library(cliqueMS)
library(SummarizedExperiment)
library(pmp)
library(heplots)
library(gplots)
library(caret)
library(ggpubr)
library(RColorBrewer)
library(oligo)
library(lumi)
library(affy)
library(limma)
library(genefilter)
library(pvca)
library(biosigner)
library(ROCR)
library(mogene21sttranscriptcluster.db)
library(AnnotationDbi)
library(biomaRt)

Metabolómica

Carga de los datos

En este ejemplo, utilizaremos una versión reducida del sacurine dataset.

pheno_dt <- phenoImport("D:/Sacurine/", header = TRUE, sep = "\t")
# We work with a reduced dataset to speed the computations
files <- list.files("D:/Jordi/Bioinformatica UOC/2n Semestre/TFM/Metabolomics/Sacurine/Processed", 
                    full.names = TRUE,pattern = "\\.mzX?ML$")
samp_names <- sub(pattern = "\\.mzX?ML$", replace = "", basename(files)) 
pd2 <- pheno_dt[which(pheno_dt$X %in% samp_names),]
str(pheno_dt)
#> 'data.frame':    240 obs. of  11 variables:
#>  $ X             : chr  "QC1_001" "Blanc04" "HU_neg_017" "HU_neg_018" ...
#>  $ sampleType    : chr  "pool" "blank" "sample" "sample" ...
#>  $ injectionOrder: int  14 16 17 18 19 20 22 23 25 26 ...
#>  $ batch         : chr  "ne1" "ne1" "ne1" "ne1" ...
#>  $ subset        : int  1 1 1 0 0 0 0 1 0 1 ...
#>  $ full          : int  1 1 1 1 1 1 1 1 1 1 ...
#>  $ sampling      : int  NA NA 1 1 1 1 1 2 2 2 ...
#>  $ osmolality    : int  NA NA 311 870 482 638 428 954 392 547 ...
#>  $ age           : num  NA NA 41 34 59 34 37 41 38 52 ...
#>  $ bmi           : num  NA NA 23 21 17.1 ...
#>  $ gender        : chr  NA NA "M" "M" ...
file_dt <- metaboImport(files[-6], phenodata = pd2, injectionvar = "X")

Visualización de los datos

Podemos visualizar el total ion count (TIC) para cada muestra ordenada por el tiempo de inyección:

ggTicQuality(file_dt,filenum = 1:50, pheno_var = 2, order = TRUE)

O filtrarlo con una variable de la phenodata y visualizar solo las muestras del control de calidad:

ggTicQuality(file_dt, pheno_var = 2, pheno_filter = "pool", order = TRUE)

Podemos ver que el TIC de los controles de calidad se va reduciendo a lo largo de las inyecciones, lo que nos podria indicar una necesidad de corregir este efecto batch.

Después podemos visualizar cromatogramas y Extracted Ion Chromatogram (EIC) (chromatograma de un intervalo de m/z pequeño) con la función ggChromPlot() según si se especifica un rango de mz o no.

ggChromPlot(file_dt, filenum = 30:35, mz = 179.056, ppm = 20, pheno_var = 1, chromtype = "sum", logscale = TRUE)

Procesado de datos

Después de visualizar los datos utilizaremos XCMS para detectar picos cromatográficos y juntarlos en features. También se utilizará CAMERA o cliqueMS para hacer la anotación de estos y se incorporará al resultado de XCMS.

proc_dt <- metaboProc(object = file_dt2, polarity = "negative", groupvar = 2,
                    peakwidth = c(10,120), noise = 0, snthresh = 20, ppm = 20,
                    expandrt = 4, binsize = 0.6, minFraction = 0.4, bw = 30,
                    annotation = "cliqueMS", cliqsamp = 1)

Para agilizar el procesado, cargaremos el resultado de la función anterior directamente:

proc_dt <- readRDS(file = "D:/Jordi/Bioinformatica UOC/2n Semestre/TFM/Metabolomics/Sacurine/Processed/proc_dt.rds")

Esta funcion produce una lista con dos o tres items dependiendo de si se ha hecho anotación. El primero es el objeto XCMSexp con todos los datos del experimento y su procesado, el segundo, el objeto de anotación y, el último, la matriz de features en formato SummarizedExperiment.

proc_dt
#> [[1]]
#> MSn experiment data ("XCMSnExp")
#> Object size in memory: 18.28 Mb
#> - - - Spectra data - - -
#>  MS level(s): 1 
#>  Number of spectra: 65591 
#>  MSn retention times: 0:1 - 19:58 minutes
#> - - - Processing information - - -
#> Data loaded [Tue Apr 06 19:23:16 2021] 
#>  MSnbase version: 2.15.7 
#> - - - Meta data  - - -
#> phenoData
#>   rowNames: 1 2 ... 121 (35 total)
#>   varLabels: X sampleType ... gender (11 total)
#>   varMetadata: labelDescription
#> Loaded from:
#>   [1] QC1_001.mzML...  [35] QC1_014.mzML
#>   Use 'fileNames(.)' to see all files.
#> protocolData: none
#> featureData
#>   featureNames: F01.S0001 F01.S0002 ... F35.S1533 (65591 total)
#>   fvarLabels: fileIdx spIdx ... spectrum (35 total)
#>   fvarMetadata: labelDescription
#> experimentData: use 'experimentData(object)'
#> - - - xcms preprocessing - - -
#> Chromatographic peak detection:
#>  method: centWave 
#>  72438 peaks identified in 35 samples.
#>  On average 2070 chromatographic peaks per sample.
#> Alignment/retention time adjustment:
#>  method: obiwarp 
#> Correspondence:
#>  method: chromatographic peak density 
#>  1360 features identified.
#>  Median mz range of features: 0.072889
#>  Median rt range of features: 106
#>  5576 filled peaks (on average 159.3143 per sample).
#> 
#> [[2]]
#> 
#> [[3]]
#> class: SummarizedExperiment 
#> dim: 1360 35 
#> metadata(5): '' '' '' '' ''
#> assays(1): raw
#> rownames(1360): FT0001 FT0002 ... FT1359 FT1360
#> rowData names(13): mzmed mzmin ... ms_level annotation
#> colnames(35): QC1_001.mzML Blanc04.mzML ... HU_neg_208.mzML
#>   QC1_014.mzML
#> colData names(11): X sampleType ... bmi gender

Batch correction

Para corregir el batch effect con QCs, utilizaremos la función batchNormalization() con el método qcnorm. Esta función utilizará la función QCRS() del paquete pmp para normalizar el efecto dentro de cada tanda y entre tandas fitando curvas polinómicas robustas.

Antes de esta corrección, es recomendable aplicar una substracción de blanco y eliminar features con alta intensidad en el blanco respecto la muestra. También eliminaremos features y muestras con alto número de valores ausentes NA.

feature_blank_na_rds <- metabFeatureFilter(proc_dt[[3]], groupvar = 2, 
                                    blankfilt = TRUE, nafilter = TRUE, 
                                    blankFoldChange = 2, blankname = "blank",
                                    naratioThr = 0.9, naratioMethod = "QC", 
                                    cvqcfilt = TRUE, cvqc_thr = 40, 
                                    qcname = "pool", sampfilter = TRUE, 
                                    maxmv = 0.4)
nrow(proc_dt[[3]]) - nrow(feature_blank_na_rds)
#> [1] 757
ncol(proc_dt[[3]]) - ncol(feature_blank_na_rds)
#> [1] 7

se han eliminado 757 y 7 muestras

library(pmp)
feature_batch <- batchNormalization(features = feature_blank_na_rds, 
                                    method = "qcnorm", injectionorder = 1:28, 
                                    groups = 2, qcname = "pool", batchnum = 4)

feature_batch <- metabFeatureFilter(feature_batch, groupvar = 2, nafilter = TRUE, 
                                    naratioThr = 0.3, naratioMethod = "across")

Ahora podemos obervar el efecto de la normalización utilizando las funciónes featurebatchQc() que nos muestra la intensidad total de las features a lo largo de las inyecciones y de PCA, que nos permite ver si se agrupan mejor los QC al normalizar:

# Feature sum of intensities plot
featurebatchQc(feature_blank_na_rds, groupvar = 2, qcname = "pool", logscale = F)

featurebatchQc(feature_batch, groupvar = 2, qcname = "pool", logscale = F)

Vemos como la normalización ha normalizado la intensidad de los QCs y ahora se mantiene más constante a lo largo de las inyecciones.

Para visualizar las diferencias en la PCA primero debemos hacer un pre-procesado de las features, se hará una normalización por la respuesta media de todos los QCs (PQN), imputacion de valores utilizando k-NN y la transformación con el logaritmo generalizado.

# PCA pre-processing and plotting
# pqn normalization, multivariate imputation and generalized logarithm are applied prior to pca plotting
feature_noblank_proc <- prePro(feature_blank_na_rds, prefuns = c("pqn", "mvImp", "glog"), 
                                method = "knn", groupvar = 2, qcname = "pool")
feature_batch_proc <- prePro(feature_batch, prefuns = c("pqn", "mvImp", "glog"), 
                                method = "knn", groupvar = 2, qcname = "pool")

raw_pca <- multipca(feature_noblank_proc, scale = FALSE)
norm_pca <- multipca(feature_batch_proc, scale = FALSE)

pcaPlot(pcdt = raw_pca, feature_noblank_proc, groupvar = 2)

pcaPlot(pcdt = norm_pca, feature_batch_proc, groupvar = 2)

Vemos como, al hacer la normalización, los controles de calidad de la PCA se clusterizan mejor.

Por último, la función sbc_plot() permite visualizar la corrección de cada feature por separado.

features_plots <- sbc_plot(feature_blank_na_rds, feature_batch, classes = colData(feature_batch)[[2]], qc_label = "pool", 
                           batch = rep(1, 28), output = NULL)
features_plots[1]
#> [[1]]

features_plots[5]
#> [[1]]

features_plots[10]
#> [[1]]

Antes de continuar con los análisis, querremos normalizar los datos originales (para no añadir los valores de la imputacion con el k-NN) y eliminar los controles de calidad.

feature_blank_na_rds <- prePro(feature_blank_na_rds, prefuns = c("pqn"), groupvar = 2, qcname = "pool")
feature_batch <- prePro(feature_batch, prefuns = c("pqn"), groupvar = 2, qcname = "pool")

feature_noblank_proc_noqc <- metabFeatureFilter(feature_blank_na_rds, groupvar = 2, 
                                                nafilter = TRUE, naratioThr = 0.4, 
                                                naratioMethod = "across", 
                                                sampfilter = TRUE, filtername = "pool")
feature_batch_no_qc <- metabFeatureFilter(feature_batch, groupvar = 2, 
                                            sampfilter = TRUE, nafilter = TRUE, naratioThr = 0.4,
                                            naratioMethod = "across", filtername = "pool")

Basic feature selection

Después del pre-procesado y la normalizacion es probable que tengamos gran cantidad de features y queramos filtrarlas para obtener las más significativas para el experimento, la función metabFeatureFilter(), que ya hemos utilizado, permite eliminar features que no corresponden a M0 (no son picos correspondientes a isotopos), no tienen una anotación o no cumplen algun criterio como, por ejemplo, tener un alto coeficiente de variación (CV).

Para decidir un criterio para el CV límite, podemos utilizar la función featureVarPlot() para ver su distribución de valores en las features:

cvFun <- function(x){sd(x, na.rm = T)/mean(x, na.rm = T)}
featureVarPlot(feature_batch_no_qc, varfun = cvFun)

Las lineas marcan los percentiles 90 y 95 de la distribución y podemos escoger entre escoger un valor absoluto límite del CV o el percentil mínimo deseado de la distribución.

selected_features <- metabFeatureFilter(feature_batch_no_qc, varfilter = TRUE, 
                                        varfun = cvFun, varquant = TRUE, 
                                        varthr = 0.7, groupvar = 2)

Feature comparison

Después de haber reducido el nombre de features a las más interesantes, las compararemos utilizando métodos univariantes. Primero, debemos decidir que test utilizaremos, el t-test lo utilizaremos cuando tengamos distribuciones normales o podamos asumir que la media es un buen estimador de las distribucion y el test de mann-whitney cuando no. Para ello, deberíamos hacer algun test de normalidad (shapiro-wilkins), QQ-plots o boxplots de las distribuciones para cada grupo. Sin embargo, como tenemos una cantidad muy pequeña de muestras, utilizaremos los test no parámetricos para este ejemplo:

group_comps_raw <- groupComp(feature_noblank_proc_noqc, groupvar = "gender", 
                            test = "wilcox", adj.method = "fdr")
group_comps_batch <- groupComp(feature_batch_no_qc, groupvar = "gender", 
                            test = "wilcox", adj.method = "fdr")
group_comps_batch_selected <- groupComp(selected_features, groupvar = "gender", 
                                        test = "wilcox", adj.method = "fdr")

Esta función genera una lista de comparaciones entre todos los factores de la variable seleccionada. Cada tabla contiene: el estadístico correspondiente, el log2 del Fold-Change, su p-valor, su p-valor ajustado por el método seleccionado y los factores comparados en la tabla (la variable en el numerador del Fold-Change corresponde a la segunda variable mencionada en esta variable).

Ahora, podemos visualizar las cmparaciones mediante volcano plots:

groupFeatureVolcano(comptable = group_comps_raw[[1]], adj.pvalue = F)

groupFeatureVolcano(comptable = group_comps_batch[[1]], adj.pvalue = F) groupFeatureVolcano(comptable = group_comps_batch_selected[[1]], adj.pvalue = F)

O filtrar features utilizando valores límite de p-value (ajustados o no) y log2(fold-change) con la función groupFeatureFilter():

selected_0_features <- groupFeatureFilter(selected_features, 
                                        comptable = group_comps_batch_selected[[1]], 
                                        pvalthr = 0.05, logFCthr = 1, 
                                        padjusted = T)
selected_0_features
#> class: SummarizedExperiment 
#> dim: 0 18 
#> metadata(7): '' '' ... original_data_structure processing_history
#> assays(1): raw
#> rownames(0):
#> rowData names(26): mzmed mzmin ... fraction fraction_flags
#> colnames(18): HU_neg_017.mzML HU_neg_028.mzML ... HU_neg_185.mzML
#>   HU_neg_204.mzML
#> colData names(14): X sampleType ... filter_samples_by_mv_flags pqn_coef
selected_features <- groupFeatureFilter(selected_features, 
                                        comptable = group_comps_batch_selected[[1]], 
                                        pvalthr = 0.05, logFCthr = 1, 
                                        padjusted = F)

selected_features
#> class: SummarizedExperiment 
#> dim: 13 18 
#> metadata(7): '' '' ... original_data_structure processing_history
#> assays(1): raw
#> rownames(13): FT0034 FT0584 ... FT1261 FT1262
#> rowData names(26): mzmed mzmin ... fraction fraction_flags
#> colnames(18): HU_neg_017.mzML HU_neg_028.mzML ... HU_neg_185.mzML
#>   HU_neg_204.mzML
#> colData names(14): X sampleType ... filter_samples_by_mv_flags pqn_coef

Vemos que, utilizando los p-valores corregidos, no encontramos ninguna feature significativa, por este motivo, en este ejemplo, utilizaremos los p-valores no corregidos aunque, los p-valores corregidos deberían ser siempre utilizados.

Para finalizar esta sección, las features serán comparadas con boxplots utilizando la función compPlot():

feature_boxplots <- compPlot(selected_features, groupvar = "gender")
feature_boxplots[[1]]

feature_boxplots[[5]]

feature_boxplots[[13]]

PCA plots

Para hacer PCA de los resultados recuperaremos las features con el efecto batch corregido, repetiremos el filtraje de las features sin eliminar los QC; las procesaremos utilizando la normalización PQN, la imputación de valores missing con un k-NN y la transformación del logaritmo generalizado (glog); eliminaremos los QCs y crearemos los plots.

feature_batch_sel <- metabFeatureFilter(feature_batch, varfilter = TRUE, 
                                        varfun = cvFun, varquant = TRUE, 
                                        varthr = 0.7, groupvar = 2, nafilter = TRUE, 
                                        naratioThr = 0.4, naratioMethod = "across")

feature_batch_sel_proc <- prePro(feature_batch_sel, 
                                prefuns = c("pqn", "mvImp", "glog"), 
                                method = "knn", groupvar = 2, qcname = "pool")
feature_batch_sel_noqc <- metabFeatureFilter(feature_batch_sel_proc, groupvar = 2, 
                                        sampfilter = TRUE, filtername = "pool")
feature_batch_sel_noqc <- groupFeatureFilter(feature_batch_sel_noqc, 
                                        comptable = group_comps_batch_selected[[1]], 
                                        pvalthr = 0.05, logFCthr = 1, 
                                        padjusted = F)
feature_sel_pca <- multipca(feature_batch_sel_noqc, groupvar = 2, scale = FALSE)
pcaPlot(feature_sel_pca, feature_batch_sel_noqc, groupvar = "gender")

Aparte del plot de la PCA, podemos querer visualizar los loadings, que nos indican la contribución de cada feature al componente correspondiente:

loadingsPlot(feature_sel_pca, pc = 1)

Para finalizar este apartado, podemos ver como se clusterizan las features en un heatmap con un dendograma en el eje de las muestras:

library(gplots)
p <- heatmapPlot(feature_batch_sel_noqc, groupvar = "gender")

Feature Selection for Machine Learning

Si volvemos depués de haber normalizado los datos por el efecto batch, nuestro objetivo puede ser el de hacer predicciones basadas en las features del dataset, por ejemplo, para predecir el sexo del paciente. Para ello, podemos utilizar el paquete biosign para seleccionar las features que son mejores predictores en tres modelos distintos: random forest, pls-da y svm. Podemos crear el modelo con la función featuresSign() y filtrar luego según el score con la función featureSelection().

feature_batch_proc <- prePro(feature_batch_no_qc, prefuns = c("mvImp"), method = "knn", groupvar = 2)

feature_sign2 <- featureSign(feature_batch_proc, groupvar = "gender")
#> No significant variable found for the selected classifier(s): 'plsda', 'randomforest', 'svm'
apply(feature_sign2@tierMC, 2, table)
#> $plsda
#> 
#>   A   B   C   E 
#>   1   1   3 589 
#> 
#> $randomforest
#> 
#>   A   E 
#>   1 593 
#> 
#> $svm
#> 
#>   E 
#> 594

Podemos ver la precisión de cada modelo y un plot con las features más significativas. En este ejemplo, solamente la feature 50 es significativa.

selected_features <- featureSelection(feature_batch_proc, feature_sign2, model = 2, scoremin = "A")

Machine Learning

Después de filtrar las features, podemos entrenar uno de los posibles modelos de machine learning:

Primero separamos el dataset en uno de validación y uno de entrenamiento: (en este ejemplo trabajaremos con las features sin filtrar por biosigner para tener una mayor cantidad de features)

set.seed(123)
train <- sample(1:18, round(ncol(feature_batch_proc)*0.7), replace = FALSE)
train_features <- feature_batch_proc[,train]
test_features <- feature_batch_proc[,-train]

Hay dos opciones disponibles para crear modelos de predicción de machine learning: utilizar los modelos creados por el paquete biosign o utilizar los modelos del paquete caret.

Para la primera opción, utilizamos la función preProcess() para hacer modificaciones en el dataset de training y obtener un objeto para pre-procesar el dataset de validación, después creamos el modelo con la función featureSign() y hacemos las predicciones con la función predict().

feature_sign2 <- featureSign(train_features, groupvar = "gender")
#> No significant variable found for the selected classifier(s): 'plsda', 'randomforest', 'svm'
apply(feature_sign2@tierMC, 2, table)
#> $plsda
#> 
#>   A   E 
#>   1 593 
#> 
#> $randomforest
#> 
#>   A   B   E 
#>   2   1 591 
#> 
#> $svm
#> 
#>   E 
#> 594
t <- table(predict(feature_sign2, newdata = t(assay(test_features)), tierMaxC = "A")[,2],
            colData(test_features)[["gender"]])
t
#>    
#>     F M
#>   F 1 3
#>   M 0 1
1 - sum(diag(t))/sum(t)
#> [1] 0.6

Para la segunda opción, la función mlFit() utilizará el modelo seleccionado del paquete caret y se obtendrá el modelo resultante, los resultados de la crossvalidation y el objeto para hacer el pre-procesado del dataset de validación:

rf_fit <- mlFit(dt = train_features, "gender", "rf", tlength = 4)
rf_fit
#> Random Forest 
#> 
#>  13 samples
#> 594 predictors
#>   2 classes: 'F', 'M' 
#> 
#> Pre-processing: centered (594), scaled (594) 
#> Resampling: Cross-Validated (10 fold, repeated 1 times) 
#> Summary of sample sizes: 12, 12, 12, 12, 12, 12, ... 
#> Resampling results across tuning parameters:
#> 
#>   mtry  Accuracy   Kappa
#>     2   0.4769231  0    
#>    13   0.5076923  0    
#>    89   0.5846154  0    
#>   594   0.5692308  0    
#> 
#> Accuracy was used to select the optimal model using  the one SE rule.
#> The final value used for the model was mtry = 89.

test_data <- predict(rf_fit$preProcess, t(assay(test_features)))

t <- table(predict(rf_fit, newdata = test_data),
        colData(test_features)[["gender"]])
t
#>    
#>     F M
#>   F 1 4
#>   M 0 0
1 - sum(diag(t))/sum(t)
#> [1] 0.8

También podemos utilizar la función mlPredictCM() para obtener la matriz de confusión correspondiente al modelo con distintos estadísticos y mlPredictROC() para obtener la curva ROC correspondiente (si el modelo permite obtener las probabilidades).

mlPredictCM(mlmod = rf_fit, newdt = test_features, 
            prepro_obj = rf_fit$preProcess, groupvar = "gender", posclass = "F")
#> Confusion Matrix and Statistics
#> 
#>        newdt_labs
#> pred_dt F M
#>       F 1 4
#>       M 0 0
#>                                           
#>                Accuracy : 0.2             
#>                  95% CI : (0.0051, 0.7164)
#>     No Information Rate : 0.8             
#>     P-Value [Acc > NIR] : 0.9997          
#>                                           
#>                   Kappa : 0               
#>                                           
#>  Mcnemar's Test P-Value : 0.1336          
#>                                           
#>             Sensitivity : 1.0             
#>             Specificity : 0.0             
#>          Pos Pred Value : 0.2             
#>          Neg Pred Value : NaN             
#>              Prevalence : 0.2             
#>          Detection Rate : 0.2             
#>    Detection Prevalence : 1.0             
#>       Balanced Accuracy : 0.5             
#>                                           
#>        'Positive' Class : F               
#> 
mlPredictROC(mlmod = rf_fit, newdt = test_features, 
            prepro_obj = rf_fit$preProcess, groupvar = "gender", posclass = "F")

Transcriptómica

Importación de los datos

direc <- "D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data"
tr_dt <- transcriImport(datapath = direc, header = TRUE, sep = ";")
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696488_WT_RT_1.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696489_WT_RT_2.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696490_WT_RT_3.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696491_KO_RT_1.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696492_KO_RT_2.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696493_KO_RT_3.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696494_WT_Cold_1.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696495_WT_Cold_2.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696496_WT_Cold_3.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696497_KO_Cold_1.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696498_KO_Cold_2.CEL
#> Reading in : D:/Jordi/Bioinformatica UOC/2n Semestre/Análisis de Datos Ómicos/StatisticalAnalysisOfMicroarrayData-master/data/GSM2696499_KO_Cold_3.CEL
tr_dt
#> GeneFeatureSet (storageMode: lockedEnvironment)
#> assayData: 1416100 features, 12 samples 
#>   element names: exprs 
#> protocolData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: exprs dates
#>   varMetadata: labelDescription channel
#> phenoData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: FileName Group ... ShortName (5 total)
#>   varMetadata: labelDescription channel
#> featureData: none
#> experimentData: use 'experimentData(object)'
#> Annotation: pd.mogene.2.1.st

Quality control

ggExprDistrPlot(tr_dt, groupvar = 1)

ggDensityPlot(tr_dt, groupvar = 1) ggPCAplot(tr_dt, groupvar = 2)

Procesado de los datos

tr_proc_dt <- procTranscript(tr_dt, annotationTable = "mogene21sttranscriptcluster.db")
#> Background correcting
#> Normalizing
#> Calculating Expression

Repetimos los plots de control de calidad después de hacer el procesado:

ggExprDistrPlot(tr_proc_dt, groupvar = 1)

ggDensityPlot(tr_proc_dt, groupvar = 1) ggPCAplot(tr_proc_dt, groupvar = 2)

Batch assesment

Vemos que porcentaje de la varianza aportan distintos factores del experimento:

featureBatchPVCA(tr_proc_dt, phenovars = 3:4, threshold = 0.3)

Como ejemplo, corregimos el efecto de la temperatura: (aunque en este caso estas diferencias nos interesan)

tr_proc_batch <- batchNormalization(tr_proc_dt, method = "covnorm", covariate = "Temperature")
featureBatchPVCA(tr_proc_batch, phenovars = 3:4, threshold = 0.3)

Vemos como, efectivamente, ahora la varianza por la temperatura es 0.

Basic feature selection

Filtramos las features con criterios básicos como un criterio de varianza (IQR) y si tienen entrada entrez

tr_proc_filt <- geneFeatureFilter(tr_proc_dt, entrez = TRUE, rem.dupEntrez = TRUE, 
                  varfilt = TRUE, var.func = IQR, varcutoff = 0.75)
tr_proc_filt
#> $eset
#> ExpressionSet (storageMode: lockedEnvironment)
#> assayData: 5991 features, 12 samples 
#>   element names: exprs 
#> protocolData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: exprs dates
#>   varMetadata: labelDescription channel
#> phenoData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: FileName Group ... ShortName (5 total)
#>   varMetadata: labelDescription channel
#> featureData: none
#> experimentData: use 'experimentData(object)'
#> Annotation: mogene21sttranscriptcluster.db 
#> 
#> $filter.log
#> $filter.log$numDupsRemoved
#> [1] 671
#> 
#> $filter.log$numLowVar
#> [1] 17973
#> 
#> $filter.log$numRemoved.ENTREZID
#> [1] 16710

Filtración con limma, anotación de tablas de comparación y comparación de grupos

Creamos las matrices de diseño y contraste y creamos un modelo de Limma con sus tablas de las comparaciones:

tr_mod_mat <- model.mat(tr_proc_filt[[1]], phenovar = 2)
tr_cont_mat <- makeContrasts(KOvsWT.COLD = KO.COLD-WT.COLD, 
                            KOvsWT.RT = KO.RT-WT.RT,
                            INT = (KO.COLD-WT.COLD) - (KO.RT-WT.RT),
                            levels = tr_mod_mat)
tr_comp_tables <- groupFeatureComp(tr_proc_filt[[1]], modelMatrix = tr_mod_mat, contrastMat = tr_cont_mat)
head(tr_comp_tables[[1]])
#>              logFC  AveExpr         t      P.Value    adj.P.Val         B
#> 17497829  2.474954 6.706683  15.71696 5.220198e-09 3.127421e-05 10.261192
#> 17407049  3.676524 4.697071  14.64453 1.117391e-08 3.347144e-05  9.706731
#> 17529307 -3.628616 3.493648 -10.95688 2.394073e-07 3.630460e-04  7.249268
#> 17373930  2.071352 5.412865  10.78024 2.832548e-07 3.630460e-04  7.105185
#> 17251565  1.991109 2.955199  10.71016 3.029928e-07 3.630460e-04  7.047237
#> 17291821 -2.434786 5.348719 -10.23434 4.835066e-07 4.153308e-04  6.641544
#>             compname
#> 17497829 KOvsWT.COLD
#> 17407049 KOvsWT.COLD
#> 17529307 KOvsWT.COLD
#> 17373930 KOvsWT.COLD
#> 17251565 KOvsWT.COLD
#> 17291821 KOvsWT.COLD
tr_ann_comp_tables <- annotateData(tableList = tr_comp_tables, anotpackage = "mogene21sttranscriptcluster.db")

Hacemos el volcano plot con 1 de las tablas de comparación:

groupFeatureVolcano(tr_ann_comp_tables[[2]])

Filtramos las features con diferencias significativas en una de las tablas:

sign_features <- groupFeatureFilter(features = tr_proc_filt[[1]], comptable = tr_ann_comp_tables[[2]], padjusted = TRUE)
sign_features
#> ExpressionSet (storageMode: lockedEnvironment)
#> assayData: 103 features, 12 samples 
#>   element names: exprs 
#> protocolData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: exprs dates
#>   varMetadata: labelDescription channel
#> phenoData
#>   rowNames: 1 2 ... 12 (12 total)
#>   varLabels: FileName Group ... ShortName (5 total)
#>   varMetadata: labelDescription channel
#> featureData: none
#> experimentData: use 'experimentData(object)'
#> Annotation: mogene21sttranscriptcluster.db

Creamos los plots con las comparaciones de cada features:

sign_features_plots <- compPlot(sign_features, groupvar = 2)
sign_features_plots[[1]]

sign_features_plots[[5]]

sign_features_plots[[10]]

Heatmap

sign_features <- annotateData(sign_features, anotpackage = "mogene21sttranscriptcluster.db")
p <- heatmapPlot(sign_features, groupvar = "Temperature")

PCA plots

tr_pca <- multipca(sign_features, scale = TRUE)
pcaPlot(tr_pca, sign_features, groupvar = 2)

loadingsPlot(tr_pca)

OmniOmics

19 de abril, 2021

Data analysis workflow examples

Metabolómica

Carga de los datos

Visualización de los datos

Procesado de datos

Batch correction

Basic feature selection

Feature comparison

PCA plots

Feature Selection for Machine Learning

Machine Learning

Transcriptómica

Importación de los datos

Quality control

Procesado de los datos

Batch assesment

Basic feature selection

Filtración con limma, anotación de tablas de comparación y comparación de grupos

Heatmap

PCA plots